<<<<<<< HEAD Case 1

Initially, we plotted the blotted weights by body weights, and we noticed that the body weights seemed to be separated into two groups.

But after color coding the chart by protocol (A,B,C,D) we noticed that the A & B protocols were grouped together and the C & D protocols were grouped together. In class, we were told that two experimental models were used, one for juvenile female rats and another for adult female ovariectomized rats; we were also told that each experimental model was associated with two protocols each. Based on this knowledge, we suspect that the gap in observed weights is likely due to A & B protocols consisting of juvenile (i.e. smaller) rats.

We then created boxplots of body and blotted weight observations across each lab. The data appeared to be right skewed, so we took the log of the response variable and replotted. The resulting boxplots appear to be more normally distributed. It also appears that all the weights vary by lab. It is also worthwhile to note that body weight and uterus weight are not measured on the same scale, but we were not given the units used for either. So, the relative size of each box plot is not to scale.

Next, we plotted the body and blotted weight observations according to the protocol that was applied. Again, we logged each of the weight values to account for data skew.

Upon plotting all combinations of dosages across groups (each group represents a different combination of dose 1 and dose 2), we notice a couple of trends. First, we see that as the quantity of dose 1 increases without dose 2, the log blotted weight increases. Secondly, when dose 2 is introduced, an increase in quantity of dose 2 leads to a decrease in log blotted weight.

group dose1 dose2 A B C D
1 0.00 0.0 89 72 54 24
2 0.00 0.0 95 72 54 24
3 0.01 0.0 84 72 54 24
4 0.03 0.0 90 72 54 23
5 0.10 0.0 96 72 54 24
6 0.30 0.0 96 72 54 24
7 1.00 0.0 96 72 54 24
8 3.00 0.0 96 72 54 24
9 10.00 0.0 96 71 54 24
10 3.00 0.1 96 72 53 24
11 3.00 1.0 96 72 54 24

Initial Model-Fitting Results

First Approach: Univariate Normal

\[ log(y) \sim \text{Norm}(\mu , \sigma^2) \]

This naive approach assumes the blotted weight follows a normal distribution with a mean centered around the mean of the log blotted weight (4.2618514) constant variance. The diagnostic plots below show that the same value is predicted for each input, the quantiles of the data do not follow a normal distrbution, and that the residuals are bimodal. We believe bimodality may be caused by the separation of rats into children and adult protocols.

Second Approach: Multivariate Normal

\[ y_i, \mu \in M_{n_i \times 1}( \Re) \\ \Sigma \in M_{n \times n}( \Re) \\ log(y_i) \sim \text{MVNorm}(\mu, \Sigma) \]

In this model,\(y_i\) is a vector of \(n_i\) observations from lab \(i\). It assumes that the log blotted weights of each lab follow approximately normal distributions with a their own means.

We provide the diagnosis plots for the first lab, Basf, and we put the others in Appendix 1.1, since R does not support plotting for multivariate regression models.

The diagnostic plots for this model indicate that the model predicts the same fitted value for each input, and that the residual quantiles do not follow a normal distribution. For the first lab, the histogram of residuals is not bimodal like they are in the simpler \(mod1\), but it does show left skew. These same problems, in addition to multimodality of the residuals, are present for the diagnostic plots for the remaining labs in the Appendix.

A weakness of this model is that, as you add more structure (blocking factors, covariates for dosage, etc.) to the model, the covariance matrix becomes more complicated, and maximum likelihood estimation becomes unwieldy. The poor fit of this model indicates that we should control for more than just the lab effect, so we consider a mixed effects model.

Third Approach: Mixed Effects

\[ y_i \in M_{n_i \times 1}( \Re) \\ \mu_i \sim \text{Norm}(\mu, \psi^2) \\ \sigma_i^2 \sim \text{IG}(\nu/2, \nu\sigma^2/2) \\ y_{ij} \sim \text{Norm}(\mu_i, \sigma^2_i) \]

Let \(y_i\) be a \(n_i \times 1\) vector of observations, where \(i\) indicates the lab and \(n_i\) is equal to the number of observations within lab \(i\). So, \(y_{ij}\) represents an individual observation for subject \(j\) within lab \(i\). We apply a Gaussian assumption here and say that each observation is normally distributed for some \(\mu_i\) and \(\sigma^2_i\).

We then add a random effect on \(\mu_i\), so \(\text{Norm}(\mu, \psi^2)\) gives us a population distribution characterizing lab-to-lab variability. This allows our model to be generalizable to future labs. Furthermore, when calculating the MLE of each \(\mu_i\), there is less variability because they are pulled in toward the grand mean \(\mu\). This introduces bias, but the decrease in varaiance results in an overall benefit for the MSE.

Analysis
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log(blot) ~ dose1 + dose2 + (1 | lab)
##    Data: dat_train
## 
##      AIC      BIC   logLik deviance df.resid 
##   4134.9   4162.9  -2062.4   4124.9     2003 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.15212 -0.85658  0.02767  0.68697  2.79594 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  lab      (Intercept) 0.1065   0.3264  
##  Residual             0.4434   0.6659  
## Number of obs: 2008, groups:  lab, 19
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  3.947951   0.077602   50.87
## dose1        0.140596   0.005203   27.02
## dose2       -0.517096   0.051522  -10.04
## 
## Correlation of Fixed Effects:
##       (Intr) dose1 
## dose1 -0.121       
## dose2 -0.051 -0.133
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log(blot) ~ dose1 + dose2 + proto + (1 | lab)
##    Data: dat_train
## 
##      AIC      BIC   logLik deviance df.resid 
##   2465.0   2509.8  -1224.5   2449.0     2000 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6148 -0.7215 -0.3028  0.7392  2.7750 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  lab      (Intercept) 0.02841  0.1686  
##  Residual             0.19331  0.4397  
## Number of obs: 2008, groups:  lab, 19
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  3.631080   0.042728   84.98
## dose1        0.142003   0.003435   41.34
## dose2       -0.519480   0.034020  -15.27
## protoB       0.037860   0.027433    1.38
## protoC       1.257924   0.030492   41.25
## protoD       1.257420   0.039567   31.78
## 
## Correlation of Fixed Effects:
##        (Intr) dose1  dose2  protoB protoC
## dose1  -0.148                            
## dose2  -0.059 -0.133                     
## protoB -0.251  0.009 -0.010              
## protoC -0.234  0.010 -0.002  0.467       
## protoD -0.159  0.009 -0.010  0.342  0.366

In our analysis, we choose to add a random effect for the lab variable, in order to account for lab-to-lab heterogeneity. Essentially, this allows us to avoid violating an independence assumption by assuming a different baseline response for each lab. First, we model these differences between individual labs by assuming random intercepts for each lab. We first include only dose1 and dose2 as fixed effects, then introduce proto as a fixed effect after confirming through anova (\(p\) < 2.2e-16, \(\chi^2\)=NA, 1675.8776) that it is a significant predictor.

Mean Fixed Effects Variance Fixed Effects
(Intercept) 3.6310802 0.0018257
dose1 0.1420033 0.0000118
dose2 -0.5194805 0.0011574
protoB 0.0378596 0.0007526
protoC 1.2579236 0.0009298
protoD 1.2574201 0.0015655
Mean Random Effects Variance Random Effects
(Intercept) 0.383359 0.0284098
Variance
Residual 0.1933111

With this full model, we see that the variability due to lab is 0.028 (in terms of variance), and variability due to non-lab sources is 0.193. And when we take a look at the fixed effects, we see that an increase in the amount of dose1 corresponds to an increase (0.142 units) in the response variable, log(blot). Furthermore, an increase in the amount of dose2 corresponds to a decrease (-0.519 units) in the response variable. Protocols B, C, and D all lead to an increase in the response variable, relative to protocol A. The diagnostic plots below show that this model fits better than the previous two. Although the residual variance decreases as the fitted values increase, the quantiles of the residuals are approximately normal.

In our previous random intercept model, we account for baseline differences between labs, but we also assume that the effects of doses is the same for each lab. After plotting the log(blot) versus dose by lab, we see that this is not a valid assumption to make since the lines are clearly not parallel. So we introduce a random slope model. Now, dose1 and dose2 can have varying slopes.

After adding the random slopes, we transformed dose 1 using a reciprocal transformation and added body weight as a predictor in order to have the best fitting model. We evaluate this model using summary statistics, diagnostic plots, and out-of-sample predictive accuracy, which we evaluate using Mean Absolute Error (\(\text{MAE} = E[|y - \hat{y}|]\)) and Root Mean Squared Error (\(\text{RMSE} = \sqrt{E[(y - \hat{y})^2]}\)). Root Mean Squared Error penalizes more for extreme error, while Mean Absolute Error simply averages all of the errors.

Mean Fixed Effects Variance Fixed Effects
(Intercept) 2.6378497 0.0728207
protoB 0.0456982 0.0003195
protoC 0.8552609 0.0094107
protoD 0.8524278 0.0100001
log(body) 0.2946495 0.0046003
Mean Random Effects Variance Random Effects
(Intercept) 3.4691702 0.9593170
I(1/(dose1 + 1/2)) 0.4969150 0.5371304
dose2 0.7531862 1.0480813
Variance
Residual 0.0797097

MAE RMSE
Random Dose + Transformations 21.79676 37.75874
Fixed Dose 35.23408 57.98964

We see that the variability due to sources other than the random effects is 0.0797 (decreased from 0.193). Also, the diagnostic plots improved, although it is notable that the residual variance decreases for larger fitted values. The fixed dose model is also much worse at predicting than the random dose model with transformations.

TODO: plot lmer results

if you plot the line for each lab and they look similar, the dose response is probably reliable if there’s a ton of variance, it’s probably not reliable because there’s a lot of variability between labs

Contributions

Nathaniel Brown made the visualizations for this report. He also organized the relevant files in a Github repository for the group to access and edit. Annie Tang compiled the group work done on EDA into a .rmd and wrote the accompanying explanations for the EDA and approaches to analysis. William Yang helped pair on EDA analysis and identify approaches to handle the data. Approaches to analysis were a joint effort by all members of the group.

Appendix

1.1 Multivariate Normal Diagnosis Plots

======= Case 1

Initially, we plotted the blotted weights by body weights, and we noticed that the body weights seemed to be separated into two groups.

But after color coding the chart by protocol (A,B,C,D) we noticed that the A & B protocols were grouped together and the C & D protocols were grouped together. In class, we were told that two experimental models were used, one for juvenile female rats and another for adult female ovariectomized rats; we were also told that each experimental model was associated with two protocols each. Based on this knowledge, we suspect that the gap in observed weights is likely due to A & B protocols consisting of juvenile (i.e. smaller) rats.

We then created boxplots of body and blotted weight observations across each lab. The data appeared to be right skewed, so we took the log of the response variable and replotted. The resulting boxplots appear to be more normally distributed. It also appears that all the weights vary by lab. It is also worthwhile to note that body weight and uterus weight are not measured on the same scale, but we were not given the units used for either. So, the relative size of each box plot is not to scale.

Next, we plotted the body and blotted weight observations according to the protocol that was applied. Again, we logged each of the weight values to account for data skew.

Upon plotting all combinations of dosages across groups (each group represents a different combination of dose 1 and dose 2), we notice a couple of trends. First, we see that as the quantity of dose 1 increases without dose 2, the log blotted weight increases. Secondly, when dose 2 is introduced, an increase in quantity of dose 2 leads to a decrease in log blotted weight.

group dose1 dose2 A B C D
1 0.00 0.0 89 72 54 24
2 0.00 0.0 95 72 54 24
3 0.01 0.0 84 72 54 24
4 0.03 0.0 90 72 54 23
5 0.10 0.0 96 72 54 24
6 0.30 0.0 96 72 54 24
7 1.00 0.0 96 72 54 24
8 3.00 0.0 96 72 54 24
9 10.00 0.0 96 71 54 24
10 3.00 0.1 96 72 53 24
11 3.00 1.0 96 72 54 24

Initial Model-Fitting Results

First Approach: Univariate Normal

\[ log(y) \sim \text{Norm}(\mu , \sigma^2) \]

This naive approach assumes the blotted weight follows a normal distribution with a mean centered around the mean of the log blotted weight (4.2618514) constant variance. The diagnostic plots below show that the same value is predicted for each input, the quantiles of the data do not follow a normal distrbution, and that the residuals are bimodal. We believe bimodality may be caused by the separation of rats into children and adult protocols.

Second Approach: Multivariate Normal

\[ y_i, \mu \in M_{n_i \times 1}( \Re) \\ \Sigma \in M_{n \times n}( \Re) \\ log(y_i) \sim \text{MVNorm}(\mu, \Sigma) \]

In this model,\(y_i\) is a vector of \(n_i\) observations from lab \(i\). It assumes that the log blotted weights of each lab follow approximately normal distributions with a their own means.

We provide the diagnosis plots for the first lab, Basf, and we put the others in Appendix 1.1, since R does not support plotting for multivariate regression models.

The diagnostic plots for this model indicate that the model predicts the same fitted value for each input, and that the residual quantiles do not follow a normal distribution. For the first lab, the histogram of residuals is not bimodal like they are in the simpler \(mod1\), but it does show left skew. These same problems, in addition to multimodality of the residuals, are present for the diagnostic plots for the remaining labs in the Appendix.

A weakness of this model is that, as you add more structure (blocking factors, covariates for dosage, etc.) to the model, the covariance matrix becomes more complicated, and maximum likelihood estimation becomes unwieldy. The poor fit of this model indicates that we should control for more than just the lab effect, so we consider a mixed effects model.

Third Approach: Mixed Effects

\[ y_{ij} \sim \beta_{0,i} + \beta_{1,i}x_{ij,d_1} + \beta_{2,i}x_{ij,d_2} + \beta_{3,i}x_{ij,p_B} + \beta_{4,i}x_{ij,p_C} + \beta_{5,i}x_{ij,p_D} + \beta_{6,i}x_{ij,log(w)} + \epsilon \\ \beta_{0:6,i} \sim N(\mu_{0:6,i}, \sigma^2_{0:6,i}) \\ \epsilon \sim N(\mu, \sigma^2) \] Let \(y_{ij}\) be the observation for log(blotted uterus weight) for subject \(x_{ij}\), the \(j\)th individual in lab \(i\). \(x_{ij,d_1}\) and \(x_{ij,d_2}\) are the values of dose1 and dose2 for subject \(x_{ij}\). \(x_{ij,p_B}\), \(x_{ij,p_C}\), and \(x_{ij,p_D}\) are dummy variables for which protocol \(x_{ij}\) was subjected to. \(x_{ij,log(w)}\) is the log(bodyweight) for \(x_{ij}\). We make the Gaussian assumption that the coefficients \(\beta\) are normally distributed according to some \(\mu_i\) and \(\sigma_i\). We add a random effect on all \(\beta_{0:5,i}\) to account for lab-to-lab variability in the intercepts and slopes on the blotted weight for the different dosages and protocols.

We start at a reduced form of this model and augment it to its full form after initial analyses.

Analysis
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log(blot) ~ dose1 + dose2 + proto + (1 | lab)
##    Data: dat_train
## 
##      AIC      BIC   logLik deviance df.resid 
##   2465.0   2509.8  -1224.5   2449.0     2000 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6148 -0.7215 -0.3028  0.7392  2.7750 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  lab      (Intercept) 0.02841  0.1686  
##  Residual             0.19331  0.4397  
## Number of obs: 2008, groups:  lab, 19
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  3.631080   0.042728   84.98
## dose1        0.142003   0.003435   41.34
## dose2       -0.519480   0.034020  -15.27
## protoB       0.037860   0.027433    1.38
## protoC       1.257924   0.030492   41.25
## protoD       1.257420   0.039567   31.78
## 
## Correlation of Fixed Effects:
##        (Intr) dose1  dose2  protoB protoC
## dose1  -0.148                            
## dose2  -0.059 -0.133                     
## protoB -0.251  0.009 -0.010              
## protoC -0.234  0.010 -0.002  0.467       
## protoD -0.159  0.009 -0.010  0.342  0.366

In our analysis, we choose to add a random effect for the lab variable, in order to account for lab-to-lab heterogeneity. Essentially, this allows us to avoid violating an independence assumption by assuming a different baseline response for each lab. First, we model these differences between individual labs by assuming random intercepts for each lab. We first include only dose1 and dose2 as fixed effects, then introduce proto as a fixed effect after confirming through anova (\(p\) < 2.2e-16, \(\chi^2\)=NA, 1675.8776) that it is a significant predictor.

Mean Fixed Effects Variance Fixed Effects
(Intercept) 3.6310802 0.0018257
dose1 0.1420033 0.0000118
dose2 -0.5194805 0.0011574
protoB 0.0378596 0.0007526
protoC 1.2579236 0.0009298
protoD 1.2574201 0.0015655
Mean Random Effects Variance Random Effects
(Intercept) 0.3833589 0.0284098
Variance
Residual 0.1933111

With this full model, we see that the variability due to lab is 0.028 (in terms of variance), and variability due to non-lab sources is 0.193. And when we take a look at the fixed effects, we see that an increase in the amount of dose1 corresponds to an increase (0.142 units) in the response variable, log(blot). Furthermore, an increase in the amount of dose2 corresponds to a decrease (-0.519 units) in the response variable. Protocols B, C, and D all lead to an increase in the response variable, relative to protocol A. The diagnostic plots below show that this model fits better than the previous two. Although the residual variance decreases as the fitted values increase, the quantiles of the residuals are approximately normal.

In our previous random intercept model, we account for baseline differences between labs, but we also assume that the effects of doses is the same for each lab. After plotting the log(blot) versus dose by lab, we see that this is not a valid assumption to make since the lines are clearly not parallel. So we introduce a random slope model. Now, dose1 and dose2 can have varying slopes.

After adding the random slopes, we transformed dose 1 using a reciprocal transformation and added body weight as a predictor in order to have the best fitting model. We evaluate this model using summary statistics, diagnostic plots, and out-of-sample predictive accuracy, which we evaluate using Mean Absolute Error (\(\text{MAE} = E[|y - \hat{y}|]\)) and Root Mean Squared Error (\(\text{RMSE} = \sqrt{E[(y - \hat{y})^2]}\)). Root Mean Squared Error penalizes more for extreme error, while Mean Absolute Error simply averages all of the errors.

Mean Fixed Effects Variance Fixed Effects
(Intercept) 2.6378477 0.0728207
protoB 0.0456982 0.0003195
protoC 0.8552608 0.0094107
protoD 0.8524277 0.0100001
log(body) 0.2946495 0.0046003
Mean Random Effects Variance Random Effects
(Intercept) 3.4691708 0.9593174
I(1/(dose1 + 1/2)) 0.4969139 0.5371285
dose2 0.7531864 1.0480777
Variance
Residual 0.0797097

MAE RMSE
Random Dose + Transformations 21.79676 37.75874
Fixed Dose 35.23408 57.98964

We see that the variability due to sources other than the random effects is 0.0797 (decreased from 0.193). Also, the diagnostic plots improved, although it is notable that the residual variance decreases for larger fitted values. The fixed dose model is also much worse at predicting than the random dose model with transformations.

TODO: plot lmer results

if you plot the line for each lab and they look similar, the dose response is probably reliable if there’s a ton of variance, it’s probably not reliable because there’s a lot of variability between labs

Contributions

Nathaniel Brown made the visualizations for this report. He also organized the relevant files in a Github repository for the group to access and edit. Annie Tang compiled the group work done on EDA into a .rmd and wrote the accompanying explanations for the EDA and approaches to analysis. William Yang helped pair on EDA analysis and identify approaches to handle the data. Approaches to analysis were a joint effort by all members of the group.

Appendix

1.1 Multivariate Normal Diagnosis Plots

>>>>>>> 8f5534e6523317996c2731e2d5755cdd381d6420